Votre projet RStudio doit être structuré de la manière suivante :
├── data
│ └── base_cc_comparateur.xls
│ └── Occitanie.gpkg
├── exo1.R
└── exercice.RprojVous devez enregistrer votre script régulièrement dans le fichier exo1.R.
Source du fichier Occitanie.gpkg : ADMIN EXPRESS COG édition 2019, IGN
Source du fichier base_cc_comparateur.xls : Base comparateur de territoires 2019, INSEE
sf et la fonction st_read().st_crs(). S’agit-il de données projetée ? Si ce n’est pas le cas, transformez la couche des communes dans la projection française (Lambert 93, EPSG : 2154) avec la fonction st_transform().library(sf)## Linking to GEOS 3.7.1, GDAL 3.1.2, PROJ 7.1.0
occ_raw <- st_read(dsn = "data/Occitanie.gpkg", stringsAsFactors = FALSE)## Reading layer `Occitanie' from data source `/home/tim/Documents/prz/carto_avec_r_exo/data/Occitanie.gpkg' using driver `GPKG'
## Simple feature collection with 4454 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -0.3271723 ymin: 42.33292 xmax: 4.845565 ymax: 45.04557
## geographic CRS: WGS 84
st_crs(occ_raw)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
occ <- st_transform(x = occ_raw, crs = 2154)st_geometry())Pour connaitre la liste de tous les noms ou code de région on peut utiliser la fonction unique().
unique(occ$INSEE_DEP)## [1] "32" "31" "82" "12" "81" "30" "11" "65" "46" "09" "34" "66" "48"
com31 <- occ[occ$INSEE_DEP == "31", ]
plot(st_geometry(com31))st_union ().aggregate() pour regrouper les polygones et calculer les sommes des populations communales.reg76 <- st_union(occ)
dep76 <- aggregate(occ[,"POPULATION"], by = list(occ$INSEE_DEP), sum)
plot(st_geometry(occ), col = "lightblue1", border = "white", lwd = .5)
plot(st_geometry(dep76), col = NA, border = "lightblue2", lwd = 1, add = TRUE)
plot(reg76, col = NA, border = "lightblue3", lwd = 2, add = TRUE)st_buffer(). Quel est le code INSEE de la commune de Toulouse?toulouse <- com31[com31$INSEE_COM == "31555", ]
toulouse_buffer <- st_buffer(toulouse, 10000)
plot(st_geometry(com31), lwd = .5)
plot(st_geometry(toulouse), col = "darkblue", add = TRUE)
plot(st_geometry(toulouse_buffer), col = NA, lwd = 3, border = "red",
add = TRUE)Déterminez quelles communes de la Haute-Garonne intersectent le buffer créé.
st_intersects().com31$in_buffer <- st_intersects(x = com31, toulouse_buffer, sparse = FALSE)
head(com31, 2)## Simple feature collection with 2 features and 12 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 504779 ymin: 6198834 xmax: 565151 ymax: 6270449
## projected CRS: RGF93 / Lambert-93
## ID STATUT INSEE_COM INSEE_ARR INSEE_DEP
## 4 COMMUNE_0000000009762857 Commune simple 31046 2 31
## 5 COMMUNE_0000000009760450 Commune simple 31547 1 31
## INSEE_REG CODE_EPCI NOM_COM_M POPULATION TYPE NOM_COM
## 4 76 200072635 BAREN 12 COM Baren
## 5 76 200068641 SEYSSES 8787 COM Seysses
## geom in_buffer
## 4 MULTIPOLYGON (((504930 6199... FALSE
## 5 MULTIPOLYGON (((556904 6267... TRUE
Affichez/superposez toutes les couches géographiques créees :
Jouez sur les styles pour les différencier…
plot(reg76, col = NA, border = "grey50", lwd = 2, bg = "lightyellow")
plot(st_geometry(com31), col = "#aec8f2", border = "white", lwd = .5,
add = TRUE)
plot(st_geometry(com31[com31$in_buffer == TRUE, ]),col = "red", lwd = .5,
border = "white", add = TRUE)
plot(st_geometry(toulouse), col = "darkblue", border = "black", lwd = 1,
add = TRUE)
plot(st_geometry(toulouse_buffer), col = NA, border = "black", lwd = 2,
lty = 2, add = TRUE)
plot(st_geometry(dep76), col = NA, border = "grey50", add = TRUE)Créer un objet sf (points) contenant la localisation de la préfécture de région à Toulouse.
st_points(), puis un objet sfc avec avec st_sfc() et finalement un objet sf avec st_sf().# Sur OSM : https://www.openstreetmap.org/search?whereami=1&query=43.59825%2C1.45047#map=19/43.59825/1.45047
# Position de la préfecture 43.59825, 1.45047
library(sf)
pref_pt <- st_point(c(1.45047, 43.59825))
pref_sfc <- st_sfc(pref_pt, crs = (4326))
pref <- st_sf(name = "Préfecture", geometry = pref_sfc)
pref## Simple feature collection with 1 feature and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 1.45047 ymin: 43.59825 xmax: 1.45047 ymax: 43.59825
## geographic CRS: WGS 84
## name geometry
## 1 Préfecture POINT (1.45047 43.59825)
Calculez une matrice de distance entre la préfecture et les centroïdes des communes du département Haute-Garonne.
Ajouter cette distance dans une nouvelle colonne dist_pref dans com31. Pour cela, utiliser les fonction st_centroid() et st_distance().
N’oubliez pas de vérifier les projections utilisées avec st_crs(), au besoin modifiez les avec st_transform().
st_crs(pref)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
st_crs(com31)## Coordinate Reference System:
## User input: EPSG:2154
## wkt:
## PROJCRS["RGF93 / Lambert-93",
## BASEGEOGCRS["RGF93",
## DATUM["Reseau Geodesique Francais 1993",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4171]],
## CONVERSION["Lambert-93",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",46.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",3,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",44,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",700000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",6600000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["France"],
## BBOX[41.15,-9.86,51.56,10.38]],
## ID["EPSG",2154]]
identical(st_crs(pref), st_crs(com31))## [1] FALSE
pref <- st_transform(pref, st_crs(com31))
identical(st_crs(pref), st_crs(com31))## [1] TRUE
# Centroides des communes
com31_centro <- st_centroid(st_geometry(com31))
com31$dist_pref <- st_distance(com31_centro, pref)
head(com31, 2)## Simple feature collection with 2 features and 13 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 504779 ymin: 6198834 xmax: 565151 ymax: 6270449
## projected CRS: RGF93 / Lambert-93
## ID STATUT INSEE_COM INSEE_ARR INSEE_DEP
## 4 COMMUNE_0000000009762857 Commune simple 31046 2 31
## 5 COMMUNE_0000000009760450 Commune simple 31547 1 31
## INSEE_REG CODE_EPCI NOM_COM_M POPULATION TYPE NOM_COM
## 4 76 200072635 BAREN 12 COM Baren
## 5 76 200068641 SEYSSES 8787 COM Seysses
## geom in_buffer dist_pref
## 4 MULTIPOLYGON (((504930 6199... FALSE 104829.70 [m]
## 5 MULTIPOLYGON (((556904 6267... TRUE 17211.48 [m]
Cartographiez distance de chaque commune du département à la préfecture avec la fonction plot().
# Le plus simple
plot(com31["dist_pref"])# Avec un peu de paramétrage
plot(com31["dist_pref"], main = "Distance à la préfecture (en mètres)",
pal = hcl.colors(13,"Turku", rev = TRUE), border = NA,
key.pos = 1,key.width = .15, key.length = .75, graticule = TRUE,
reset = FALSE)
plot(st_geometry(pref), pch = 20, col = "red", add = TRUE )osmdata extrayez l’ensemble des restaurants présents dans le département.com46 <- occ[occ$INSEE_DEP == "46", ]
plot(st_geometry(com46))library(osmdata)## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
# Définition d'une bounding box
q <- opq(bbox=st_bbox(st_transform(com46,4326)))
# Extraction des restaurants
res <- add_osm_feature(opq = q, key = 'amenity', value = "restaurant")
res.sf <- osmdata_sf(res)
res.sf.pts <- res.sf$osm_points[!is.na(res.sf$osm_points$amenity),]
resto <- st_transform(res.sf.pts, st_crs(com46))
resto <- st_intersection(resto, st_geometry(com46))## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# Affichage des restaurants
plot(st_geometry(com46), col="darkseagreen3", border="darkseagreen4")
plot(st_geometry(resto), add=TRUE, pch=20, col = "#330A5FFF", cex = 0.5)# Compter les restaurants par communes
inter <- st_intersects(x = com46, y = resto)
com46$nresto <- sapply(inter, length)
# communes sans restos
com46noresto <- com46[com46$nresto==0, ]
plot(st_geometry(com46), col="darkseagreen3", border="darkseagreen4")
plot(st_geometry(com46noresto), col = 'red', add = TRUE)
plot(st_geometry(resto), add=TRUE, pch=20, col = "blue", cex = 1)# index du restaurant le plus proche
index <- st_nearest_feature(x = st_centroid(com46noresto),
y = resto)## Warning in st_centroid.sf(com46noresto): st_centroid assumes attributes are
## constant over geometries of x
# distance au plus proche
com46noresto$dresto <- st_distance(x = st_centroid(com46noresto),
y = resto[index, ],
by_element = TRUE)## Warning in st_centroid.sf(com46noresto): st_centroid assumes attributes are
## constant over geometries of x
# Affichage de la carte
plot(com46noresto['dresto'], reset = F,
main = "Distance au restaurant le plus proche")
plot(st_geometry(com46), col=NA, add= TRUE)
plot(st_geometry(resto), add=TRUE, pch=20, col = "red", cex = 1)Dans le même projet RStudio créez un deuxième script nommé exo2.R.
Votre projet RStudio doit maintenant être structuré de la manière suivante :
├── data
│ └── base_cc_comparateur.xls
│ └── Occitanie.gpkg
├── exo1.R
├── exo2.R
└── exercice.Rprojsf et la fonction st_read() pour importer les données.st_crs().readxl et la fonction read_excel() pour ouvrir la table de données correctement. Importer la table dans un objet nommé occ_df.library(sf)
occ_raw <- st_read(dsn = "data/Occitanie.gpkg", stringsAsFactors = FALSE)## Reading layer `Occitanie' from data source `/home/tim/Documents/prz/carto_avec_r_exo/data/Occitanie.gpkg' using driver `GPKG'
## Simple feature collection with 4454 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -0.3271723 ymin: 42.33292 xmax: 4.845565 ymax: 45.04557
## geographic CRS: WGS 84
st_crs(occ_raw)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
occ <- st_transform(x = occ_raw, crs = 2154)
library(readxl)
occ_df <- read_excel(path = "data/base_cc_comparateur.xls", sheet = 1, skip = 5)
head(occ_df, 2)## # A tibble: 2 x 36
## CODGEO LIBGEO REG DEP P16_POP P11_POP SUPERF NAIS1116 DECE1116 P16_MEN
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 01001 L'Abe… 84 01 767 780 16.0 41 25 306
## 2 01002 L'Abe… 84 01 243 234 9.15 21 7 101
## # … with 26 more variables: NAISD18 <dbl>, DECESD18 <dbl>, P16_LOG <dbl>,
## # P16_RP <dbl>, P16_RSECOCC <dbl>, P16_LOGVAC <dbl>, P16_RP_PROP <dbl>,
## # NBMENFISC16 <dbl>, PIMP16 <dbl>, MED16 <dbl>, TP6016 <dbl>,
## # P16_EMPLT <dbl>, P16_EMPLT_SAL <dbl>, P11_EMPLT <dbl>, P16_POP1564 <dbl>,
## # P16_CHOM1564 <dbl>, P16_ACT1564 <dbl>, ETTOT15 <dbl>, ETAZ15 <dbl>,
## # ETBE15 <dbl>, ETFZ15 <dbl>, ETGU15 <dbl>, ETGZ15 <dbl>, ETOQ15 <dbl>,
## # ETTEF115 <dbl>, ETTEFP1015 <dbl>
Créez un nouvel objet sf à partir d’une séléction par attribut`:
plot(st_geometry(...)).Pour connaitre la liste de tous les noms ou code de région voir la fonction unique()
unique(occ$INSEE_DEP)## [1] "32" "31" "82" "12" "81" "30" "11" "65" "46" "09" "34" "66" "48"
com31 <- occ[occ$INSEE_DEP=="31", ]
plot(st_geometry(com31))merge(). Quel est l’identifiant commun entre la table géo et la table attributaire?Réalisez huit cartes thématiques finalisées.
Chacune des cartes devra comprendre les éléments d’habillage et de mise en page nécessaires leur compréhension (titre, légende, sources…)
Parmi ces huits cartes deux échelles différentes (échelle de la région, échelle d’un département) et deux découpages territoriaux différents (les communes, les départements) doivent être utilisés.
Vous devez réaliser les huit types de cartes suivants :
Exemple :
Nous vous conseillons d’utiliser le package cartography pour l’ensemble de vos réalisations. Néanmoins l’utilisation d’autres packages comme ggplot2 ou tmap est possible.
Vous devrez nous rendre pour le xx février 2020 un script permettant de construire ces huit cartes dans un fichier dont le nom est constitué de votre prénom et de votre nom (Prenom_Nom.R).
Votre script doit pouvoir fonctionner une fois inséré dans un projet RStudio structuré de la manière suivante :
├── data
│ ├── base_cc_comparateur.xls
│ └── Occitanie.gpkg
├── Prenom_Nom.R
└── exercice.RprojVous serez évalués de la manière suivante :
Pour répondre à cet exercice vous devrez créer au moins deux nouvelles variables. Une variable quantitative relative :
Une variable qualitative (n’utilisez par cette variables dans votre exercice ) :